\input{euler.tex}

\newcommand{\vA}{\mathbf{A}}
\newcommand{\vB}{\mathbf{B}}
\newcommand{\vI}{\mathbf{I}}
\newcommand{\vS}{\mathcal{S}}
\newcommand{\vT}{\mathbf{T}}
\newcommand{\vone}{\mathbf{1}}

\begin{document}

\problem[368]{A Kempner-Like Series}

It is well known that the harmonic series is divergent. That is,
\[
\frac11 + \frac12 + \frac13 + \frac14 + \cdots = \infty .
\]

However, if we omit from this series every term where the denominator contains the digit 9, the series remarkably enough converges to approximately 22.9206766193. This modified harmonic series is called the \emph{Kempner series}.

Let us now consider another modified harmonic series by omitting from the harmonic series every term where the denominator contains three or more consecutive equal digits. One can verify that out of the first 1200 terms of the harmonic series, only 20 terms will be omitted. These 20 omitted terms are:
\[
\frac{1}{111}, \cdots, \frac{1}{999}, \frac{1}{1000}, \frac{1}{1110}, \frac{1}{1111}, \cdots, 
\frac{1}{1119} .
\]
This series converges as well.

Find the value the series converges to. Give your answer rounded to 10 digits behind the decimal point.

\solution

Schmelzer and Baillie (2008) published an algorithm that computes the convergence value of the generalized Kempner series efficiently. Below is a summary of some key points in their method, followed by an application to this particular problem.

Let $\Psi$ denote the sum of the generalized Kempner sequence. Let $s$ denote the denominator of a term in the sequence, and let $\vS$ denote the set of all the denominators. Then our goal is to find the value of
\begin{equation}
\Psi = \sum_{s \in \vS} \frac{1}{s}. \label{eq:368.1}
\end{equation}
Since the Kempner series converges extremely slowly, a direct computation of the sum using \eqref{eq:368.1} is infeasible. An efficient algorithm makes use of several key techniques to achieve fast convergence.

First, group the terms by the number of digits in their denominators, and sum each group in order. Let $\vS_i$ denote the set of denominators with $i$ digits, and let $\Psi_i$ be the sum of those terms. That is,
\begin{equation}
\Psi_i = \sum_{s \in \vS_i} \frac{1}{s}. \label{eq:368.2}
\end{equation}
It turns out that while the number of terms in each group grows exponentially, the sum of each group actually decreases geometrically when $i$ is large.

Second, partition the terms into $n$ classes, $\vS_i^j (1 \le j \le n)$, and define an $n$-by-10 transformation matrix $\vT$, such that appending the digit $d$ to an arbitrary element in class $\vS_i^j$ results in an element in class $\vS_{i+1}^{T(j,d)}$, regardless of $i$. For example, in this problem, we could partition the terms into 20 classes, namely
\begin{center}
\begin{tabular}{ c l }
$\vS^1$    & terms that end with 0 but not 00 \\
           & $\cdots$ \\
$\vS^{10}$ & terms that end with 9 but not 99 \\
$\vS^{11}$ & terms that end with 00 \\
           & $\cdots$ \\
$\vS^{20}$ & terms that end with 99 \\
\end{tabular}
\end{center}
And define the transformation matrix
\[
\vT = .
\]
In particular, if appending digit $d$ to class $j$ leads to a term that is omitted, then define $T(j,d) = 0$. With this setup, the sum can be written as
\[
\Psi = \sum_{i=1}^{\infty} \sum_{j=1}^n \sum_{s \in \vS_i^j} \frac{1}{s} .
\]

At this point it's useful to define $\Psi_i^j(k) \equiv \sum_{s \in \vS_i^j} s^{-k}$. It is also useful to define the indicator function $f_{jlm} \equiv \vone_{T(l,m)=j}$ to indicate that an element class in $j$ may be obtained by appending the digit $m$ to an element in class $l$. Then we have the following recurrence relation
\begin{align}
\Psi_{i}^j(k) \equiv \sum_{s \in \vS_i^j} s^{-k}
&= \sum_{T(l,m)=j} \sum_{s \in \vS_{i-1}^l} (10 s + m)^{-k} , \notag \\
&= \sum_{l=1}^n \sum_{m=0}^9 f_{jlm} \sum_{s \in \vS_{i-1}^l} (10 s + m)^{-k} . \notag
\end{align}

Next, we make the critical step by expanding $(10s+m)^{-k}$ using negative bonimal series
\[
(10s+m)^{-k} = \sum_{w=0}^{\infty} (-1)^w \binom{k+w-1}{w} (10s)^{-k-w} m^w .
\]
Then we have
\begin{equation}
\Psi_i^j(k) = \frac{1}{10^k} \sum_{l=1}^n \sum_{m=0}^9 f_{jlm} \sum_{w=0}^{\infty} \left(-\frac{m}{10}\right)^w \binom{k+w-1}{w} \Psi_{i-1}^l (k+w) . \label{eq:368.20}
\end{equation}
And the required sum is
\begin{equation}
\Psi = \sum_{i=1}^{\infty} \sum_{j=1}^n \Psi_i^j(1). \label{eq:368.21}
\end{equation}

Up till now, equation \eqref{eq:368.21} is still exact. However, it involves two infinite sums: one in $w$ and one in $i$. Therefore, to compute it, we need to truncate small terms in $w$ and $i$. To truncate $w$, we give a priori an integer $W$ and ignore all $\Psi_{i-1}^j(k+w)$ where $k+w > W$. To truncate $i$, ignore for $i \ge K$ all $\Psi_{i-1}^j(k+w)$ where $k+w > 1$, which simplifies \eqref{eq:368.20} as
\begin{equation}
\Psi_i^j \approx \frac{1}{10} \sum_{l=1}^n \sum_{m=0}^9 f_{jlm} \Psi_{i-1}^l . 
\label{eq:368.25}
\end{equation}
Equation \eqref{eq:368.25} can be written in matrix form as
\begin{equation}
\begin{pmatrix}
\Psi_i^1 \\ \vdots \\ \Psi_i^n
\end{pmatrix}
\approx \left[ \frac{1}{10} \sum_{m=0}^9 
\begin{pmatrix}
f_{11m} & \cdots & f_{1nm} \\
\vdots  & & \vdots \\
f_{n1m} & \cdots & f_{nnm} \\
\end{pmatrix}
\right] 
\begin{pmatrix}
\Psi_{i-1}^1 \\ \vdots \\ \Psi_{i-1}^n
\end{pmatrix}
. \label{eq:368.26}
\end{equation}
Denoting the matrix in brackets by $\vA$, the remaining terms in \eqref{eq:368.21}, starting from $i = K$, can be approximated by
\[
\sum_{i=0}^{\infty}
\begin{pmatrix}
\Psi_{K+i}^1 \\ \vdots \\ \Psi_{K+i}^n
\end{pmatrix}
\approx \sum_{i=0}^{\infty} \vA^i 
\begin{pmatrix}
\Psi_{K}^1 \\ \vdots \\ \Psi_{K}^n
\end{pmatrix}
.
\]
Note that $\sum_{i=0}^{\infty} \vA^i = (\vI - \vA)^{-1}$. Hence the sum is approximated as
\begin{equation}
\Psi \approx \sum_{i=1}^{K-1} \sum_{j=1}^n \Psi_i^j(1) + \vone'_n (\vI - \vA)^{-1} \begin{pmatrix}
\Psi_{K}^1 \\ \vdots \\ \Psi_{K}^n
\end{pmatrix} 
. \label{eq:368.30}
\end{equation}

In summary, we compute $\Psi_1(1)$ and $\Psi_2(k)$ for $1 \le k \le W$ directly using \eqref{eq:368.2}. Then we compute \eqref{eq:368.20} iteratively. After each iteration, we compute $\Psi$ using \eqref{eq:368.30}. Continue until the difference between two successive computation of $\Psi$ is smaller than a given tolerance.

\complexity

Time complexity: $\BigO(1)$.

Space complexity: $\BigO(1)$.

\answer

\reference

http://en.wikipedia.org/wiki/Kempner\_series

http://eprints.maths.ox.ac.uk/1106/1/NA-06-17.pdf

http://mathworld.wolfram.com/NegativeBinomialSeries.html

\end{document} 